Calculate functional diversity metrics

knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE, 
                      fig.height = 4, fig.width = 7)

library(raster)
library(sf)
library(tidyverse)
library(here)
library(oharac)
source(here('common_fxns.R'))

Summary

Using functional entities defined in a prior script based on a suite of categorical traits, in combination with spatial distribution of species from AquaMaps and IUCN, calculate functional diversity metrics based on Mouillot et al (2014) and code from Sebastien Villeger: functional richness, functional vulnerability, functional redundancy, and functional over-redundancy.

Methods

Summarize functional entity membership per cell across IUCN and AquaMaps

Included species are:

• coded for functional entity membership
• mapped in one (or both) of the species distribution mapsets
• not necessarily part of the trait-based vulnerability set

Identify species included in Functional Entity calculations

spp_fe <- read_csv(here('_output/func_entities/fe_species.csv'), show_col_types = FALSE) %>%
  select(species) %>%
  mutate(fe = TRUE) %>%
  distinct()
spp_am <- get_am_spp_info()  %>%
  filter(occur_cells >= 10) %>%
  select(species = sciname) %>%
  mutate(am_mapped = TRUE) %>%
  distinct()
spp_iucn <- read_csv(here('_data/iucn_spp/iucn_to_worms_match.csv'), show_col_types = FALSE) %>%
  rename(species = worms_name, iucn_mapped = mapped)

spp_vuln <- get_spp_vuln() %>%
  select(species) %>%
  mutate(vuln = TRUE) %>%
  distinct()

spp_worms <- assemble_worms() %>%
  select(species) %>%
  mutate(worms = TRUE) %>%
  distinct()

all_spp <- spp_worms %>%
  full_join(spp_vuln, by = 'species') %>%
  full_join(spp_am,   by = 'species') %>%
  full_join(spp_iucn, by = 'species') %>%
  full_join(spp_fe,   by = 'species')

spp_for_fe_metrics <- all_spp %>%
  ### names in WoRMS
  filter(worms) %>%
  ### only those spp in vulnerability project?
  # filter(vuln) %>%
  ### mapped in AquaMaps and/or IUCN
  filter(am_mapped | iucn_mapped) %>%
  ### functional traits available
  filter(fe)

# worms_mismatch <- all_spp %>% 
#   filter(is.na(worms)) %>%
#   filter()
# 
# ### all AquaMaps spp
# am_check <- get_am_spp_info() %>%
#   filter(sciname %in% worms_mismatch$species)
# table(am_check %>% select(kingdom, phylum))
# table(am_check %>% filter(kingdom == 'animalia') %>% select(class, kingdom))

Load species ranges and summarize to functional entities per cell

Here we use the IUCN species range maps rasterized to 10 km Mollweide, and AquaMaps species distributions reprojected to 10 km Mollweide. Priority: IUCN species maps, then AquaMaps with occur_cells ≥ 10.

Read in functional entity assignments and join to spatial ranges. Per cell per FE, calculate the number of species (separately for IUCN species and AquaMaps species) for later calculation into functional richness, redundancy, vulnerability.

iucn_fe_summary_file <- here_anx('func_entities', 'fe_cell_summary_iucn.csv')
# unlink(iucn_fe_summary_file)
if(!file.exists(iucn_fe_summary_file)) {
  
  spp_fe <- read_csv(here('_output/func_entities/fe_species.csv'), show_col_types = FALSE)
  
  iucn_to_map <- spp_for_fe_metrics %>%
    left_join(spp_fe, by = 'species') %>%
    filter(iucn_mapped) %>%
    select(fe_id, iucn_sid) %>%
    distinct() %>%
    arrange(iucn_sid)
  
  iucn_map_fstem <- here_anx('spp_maps_mol', 'iucn_spp_mol_%s.csv')
  iucn_map_fs <- sprintf(iucn_map_fstem, iucn_to_map$iucn_sid)

  message('Assembling IUCN species maps and binding FEs...')
  iucn_maps <- parallel::mclapply(iucn_map_fs, mc.cores = 40, 
                  FUN = function(f) {
                    data.table::fread(f) %>%
                      filter(presence != 5) %>% 
                      select(-presence)
                  }) %>%
    setNames(iucn_to_map$fe_id) %>%
    data.table::rbindlist(idcol = 'fe_id') %>%
    mutate(fe_id = as.integer(fe_id))
  
  iucn_fes <- iucn_to_map$fe_id %>% unique() %>% sort()
  
  ### summarize was having a hard time with this massive df, so
  ### try breaking it down by FE and parallel processing it
  message('Summarizing IUCN Functional Richness map...')
  iucn_fe_summary <- parallel::mclapply(iucn_fes, mc.cores = 4,
                        FUN = function(fe) { ### fe <- iucn_fes[1]
                          df <- iucn_maps %>%
                            filter(fe_id == fe)
                          message('Processing fe = ', fe, '... ', nrow(df), ' species-cell combos')
                          df_out <- df %>%
                            group_by(cell_id) %>%
                            summarize(n_spp_fe = n(),
                                      fe_id = fe,
                                      .groups = 'drop')
                          return(df_out)
                        }) %>%
    data.table::rbindlist() %>%
    mutate(src = 'iucn')
  
  rm('iucn_maps')
  
  write_csv(iucn_fe_summary, iucn_fe_summary_file)
}

am_fe_summary_file <- here_anx('func_entities', 'fe_cell_summary_am.csv')
# unlink(am_fe_summary_file)
if(!file.exists(am_fe_summary_file)) {
  
  spp_fe <- read_csv(here('_output/func_entities/fe_species.csv'), show_col_types = FALSE)
  
  am_to_map <- spp_for_fe_metrics %>%
    left_join(spp_fe, by = 'species') %>%
    filter(am_mapped) %>%
    filter(!iucn_mapped | is.na(iucn_mapped)) %>%
    select(fe_id, species) %>%
    distinct() %>%
    arrange(species)
  
  am_map_fstem <- here_anx('spp_maps_mol', 'am_spp_mol_%s.csv')
  am_map_fs <- sprintf(am_map_fstem, str_replace_all(am_to_map$species, ' ', '_'))
  
  message('Assembling Aquamaps species maps and binding FEs...')
  am_maps <- parallel::mclapply(am_map_fs, mc.cores = 40, 
                  FUN = function(f) {
                    data.table::fread(f) %>%
                      filter(prob >= 0.5) %>% 
                      select(-prob)
                  }) %>%
    setNames(am_to_map$fe_id) %>%
    data.table::rbindlist(idcol = 'fe_id') %>%
    mutate(fe_id = as.integer(fe_id))
  
  am_fes <- am_to_map$fe_id %>% unique() %>% sort()
  
  ### summarize was having a hard time with this massive df, so
  ### try breaking it down by FE and parallel processing it
  message('Summarizing Aquamaps Functional Richness map...')

  # am_fe_summary <- parallel::mclapply(am_fes, mc.cores = 2,
  am_fe_summary <- lapply(am_fes, ### for memory issues, just use a single core... ugh
                        FUN = function(fe) { ### fe <- am_fes[1]
                          df <- am_maps %>%
                            filter(fe_id == fe)
                          message('Processing fe = ', fe, '... ', nrow(df), ' species-cell combos')
                          df_out <- df %>%
                            group_by(cell_id) %>%
                            summarize(n_spp_fe = n(),
                                      fe_id = fe,
                                      .groups = 'drop')
                          return(df_out)
                        }) %>%
    data.table::rbindlist() %>%
    mutate(src = 'am')
  
  rm('am_maps')
  
  write_csv(am_fe_summary, am_fe_summary_file)
}

tot_fe_summary_file <- here_anx('func_entities/fe_cell_summary_total.csv')
### unlink(tot_fe_summary_file)
if(!file.exists(tot_fe_summary_file)) {
  am_fe_summary   <- data.table::fread(am_fe_summary_file)
  iucn_fe_summary <- data.table::fread(iucn_fe_summary_file)
  
  fe_sum_maps <- data.table::rbindlist(list(am_fe_summary, iucn_fe_summary))
  
  tot_fes <- fe_sum_maps$fe_id %>% unique() %>% sort()

  tot_fe_summary <- parallel::mclapply(tot_fes, mc.cores = 12,
                      FUN = function(fe) { ### fe <- tot_fes[1]
                        message('Summarizing across AM and IUCN for FE #', fe)
                        df <- fe_sum_maps %>%
                          filter(fe_id == fe) %>%
                          group_by(cell_id) %>%
                          summarize(n_spp_fe = sum(n_spp_fe, na.rm = TRUE),
                                    fe_id = fe,
                                    .groups = 'drop')
                        return(df)
                      }) %>%
      data.table::rbindlist()
  
  write_csv(tot_fe_summary, tot_fe_summary_file)
}

Calculate functional diversity metrics

Calculate functional diversity metrics from IUCN and AquaMaps summaries of functional entity membership per cell.

This set of calculations replicates the results of Villeger’s function but without resorting to a presence/absence matrix, instead just relying on the original AquaMaps species-cell table.

• spp_fe_df joined to am_spp_cell;
• group by cell and functional group;
• calculate number of spp per functional group;
• group by cell;
• calculate functional entity metrics -
  • Species richness
  • Functional entity richness
  • Functional vulnerability (% of FEs represented by a single spp)
  • Weighted functional vulnerability (new: vuln of an FE is \(\frac{1}{2^{n-1}}\), take average over all FE)
  • Functional redundancy (mean spp count per FE)
  • Functional over-redundancy (sum of spp per FE above mean, divided by total number of spp)

fe_metrics_file <- here_anx('func_entities/fe_metrics_per_cell.csv')
### unlink(fe_metrics_file)
if(!file.exists(fe_metrics_file)) {
  # system.time({
  fe_sum_total <- data.table::fread(tot_fe_summary_file)

  ### there are 6.5 million cells.  Chop into 500k instances and mclapply to summarize
  # cell_id_vec <- 1:ncell(raster::raster(here('_spatial/ocean_area_mol.tif')))
  chunk_size <- 500000
  n_chunks <- ceiling(ncell(raster::raster(here('_spatial/ocean_area_mol.tif'))) / chunk_size)
  # system.time({
    message('Calculating species richness map...')
    nspp_df <- parallel::mclapply(1:n_chunks, mc.cores = n_chunks,
                  FUN = function(n) { ### n <- 1
                    cell_id_min <- (n - 1) * chunk_size + 1
                    cell_id_max <- n * chunk_size
                    nspp_sum <- fe_sum_total %>%
                      filter(between(cell_id, cell_id_min, cell_id_max)) %>%
                      group_by(cell_id) %>%
                      summarize(n_spp = sum(n_spp_fe, na.rm = TRUE),
                                n_fe  = n_distinct(fe_id),
                                .groups = 'drop')
                    }) %>%
      data.table::rbindlist()
  # }) ### 30 sec
  
  # system.time({
    message('Calculating functional vulnerability map...')
    f_vuln_df <- parallel::mclapply(1:n_chunks, mc.cores = n_chunks,
                FUN = function(n) { ### n <- 1
                  cell_id_min <- (n - 1) * chunk_size + 1
                  cell_id_max <- n * chunk_size
                  f_vuln <- fe_sum_total %>%
                    filter(between(cell_id, cell_id_min, cell_id_max)) %>%
                    oharac::dt_join(nspp_df, by = 'cell_id', type = 'left') %>%
                    group_by(cell_id) %>%
                    summarize(n_fe = first(n_fe), n_spp = first(n_spp),
                              f_vuln  = sum(n_spp_fe == 1) / n_fe, 
                              f_wvuln = mean(0.5^(n_spp_fe - 1)),
                              f_red   = mean(n_spp_fe),
                              .groups = 'drop')
                }) %>%
      data.table::rbindlist()
  # }) ### 30 sec

  # system.time({
    message('Calculating functional overredundancy map...')
    f_overred_df <- parallel::mclapply(1:n_chunks, mc.cores = n_chunks,
                FUN = function(n) { ### n <- 1
                  cell_id_min <- (n - 1) * chunk_size + 1
                  cell_id_max <- n * chunk_size
                  f_overred <- fe_sum_total %>%
                    filter(between(cell_id, cell_id_min, cell_id_max)) %>%
                    oharac::dt_join(f_vuln_df, by = 'cell_id', type = 'left') %>%
                    mutate(f_over = ifelse(n_spp_fe > f_red, n_spp_fe - f_red, 0)) %>%
                    group_by(cell_id) %>%
                    summarize(f_overred = sum(f_over) / sum(n_spp_fe),
                              .groups = 'drop')
                }) %>%
      data.table::rbindlist()
    # }) ### 7 sec
    
    message('Combining functional diversity metrics maps...')

    fe_metrics_df <- f_vuln_df %>%
      oharac::dt_join(f_overred_df, by = 'cell_id', type = 'full') %>%
      mutate(across(where(is.numeric), .fns = ~round(.x, 6)))
      ### rounding to 6 dec places saves a lot of storage for tifs
  # })
    # ggplot(fe_metrics_df, aes(x = f_vuln, y = f_wvuln)) +
    #   geom_point(alpha = .3) +
    #   geom_abline(slope = 1, intercept = 0)
  
  write_csv(fe_metrics_df, fe_metrics_file)
}

Map functional diversity metrics

Create and save rasters of functional diversity metrics

fe_metrics_df <- data.table::fread(fe_metrics_file)

n_spp_rast   <- map_to_mol(fe_metrics_df, which = 'n_spp')
n_fe_rast    <- map_to_mol(fe_metrics_df, which = 'n_fe')
f_vuln_rast  <- map_to_mol(fe_metrics_df, which = 'f_vuln')
f_wvuln_rast <- map_to_mol(fe_metrics_df, which = 'f_wvuln')
f_red_rast   <- map_to_mol(fe_metrics_df, which = 'f_red')
f_overred_rast <- map_to_mol(fe_metrics_df, which = 'f_overred')

writeRaster(n_spp_rast,   here('_output/func_entities/n_spp.tif'),   
            overwrite = TRUE)
writeRaster(n_fe_rast,    here('_output/func_entities/n_fe.tif'),    
            overwrite = TRUE)
writeRaster(f_vuln_rast,  here('_output/func_entities/f_vuln.tif'), 
            overwrite = TRUE)
writeRaster(f_wvuln_rast, here('_output/func_entities/f_wvuln.tif'), 
            overwrite = TRUE)
writeRaster(f_red_rast,   here('_output/func_entities/f_red.tif'),   
            overwrite = TRUE)
writeRaster(f_overred_rast, here('_output/func_entities/f_overred.tif'), 
            overwrite = TRUE)

Examine distributions of each raster.

Plot each raster.